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We realize and study the ionic Hubbard model using an interacting two-component gas of fermionic 
atoms loaded into an optical lattice. The bipartite lattice has honeycomb geometry with a staggered 
energy-offset that explicitly breaks the inversion symmetry. Distinct density-ordered phases are 
identified using noise correlation measurements of the atomic momentum distribution. For weak 
interactions the geometry induces a charge density wave. For strong repulsive interactions we detect 
a strong suppression of doubly occupied sites, as expected for a Mott insulating state, and the 
externally broken inversion symmetry is not visible anymore in the density distribution. The local 
density distributions in different configurations are characterized by measuring the number of doubly 
occupied lattice sites as a function of interaction and energy-offset. We further probe the excitations 
of the system using direction dependent modulation spectroscopy and discover a complex spectrum, 
which we compare with a theoretical model. 
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Changes in the fundamental properties of interacting 
many-body systems are often determined by the com¬ 
petition between different energy scales, which may in¬ 
duce phase transitions. A particularly intriguing situa¬ 
tion arises when the geometry of a system sets an energy 
scale that competes with the scale given by the inter¬ 
action of its constituents. The importance of geometry 
is apparent in reduced dimensions which influences the 
interacting many-body system in its evolution from one 
phase to another [1]. A tractable approach to generic 
questions is provided by the ionic Hubbard model, which 
captures key aspects of the physics of a competing geom¬ 
etry and interactions in the charge sector. The Hamilto¬ 
nian has a staggered energy-offset on a bipartite lattice, 
such that geometry supports a band insulating charge 
density wave (CDW). Conversely, strong repulsive on¬ 
site interactions favour a Mott insulating state (MI) at 
half-filling, which does not reflect the broken symmetry 
of the underlying lattice. The model was introduced in 
the context of charge-transfer organic salts [3, 4] and has 
been proposed to explain strong electron correlations in 
ferroelectric perovskite materials [5]. Ultracold atoms 
in optical lattices are an excellent platform for studying 
competing energy scales, as they allow for tuning various 
parameters and the geometry of the Hamiltonian [6-18]. 
Here we explore the ionic Hubbard Model using ultra¬ 
cold fermions loaded into a tunable optical honeycomb 
potential. 

The ionic Hubbard model has been studied theoreti¬ 
cally in ID chains [19-24] and on the 2D square lattice 
[25-28] . More recently, these studies have been extended 
to a honeycomb lattice, motivated by possible connec¬ 
tions to superconductivity in layered nitrides [29] and 
strongly correlated topological phases [30]. We consider 


the ionic Hubbard model on a honeycomb lattice: 

H = -t Yi 4A<t + U + A E W 

(ij),cr i i£A,cr 

where cj a and c ia are the creation and annihilation op¬ 
erators of one fermion with spin a = t,| on site i and 
nia = c| a Ci a . The system is characterized by three en¬ 
ergies: the kinetic energy denoted by the tunnelling am¬ 
plitude t and summed over nearest neighbours ( ij ), the 
on-site interaction U and the staggered energy-offset be¬ 
tween sites of the A and B sub-lattice A, with A > 0. 
In addition there is a harmonic confinement in all three 
directions. All parameters of the Hamiltonian are com¬ 
puted using Wannier functions [31]. 

The interplay between the interaction energy U, the 
energy-offset A and tunnelling t leads to quantum phases 
which differ by their density ordering. The two limiting 
cases can be qualitatively understood in the atomic limit 
at half-filling. For V A the system is described by a 
MI state. For a large energy-offset A >> U, we expect a 
band insulator with staggered density and two fermions 
on lattice site B [25]. The resulting CDW pattern re¬ 
flects the broken inversion symmetry of the underlying 
geometry. We can characterize the transition by an or¬ 
der parameter Na — A^b, which is zero in the MI state 
or when A = 0, with Na(b) the total number of atoms 
on sub-lattice A(B). Fig. la provides a schematic view 
of the different scenarios. 

In order to realize the ionic Hubbard model we cre¬ 
ate a quantum degenerate cloud of 40 K as described in 
previous work [31] and detailed in [32]. We prepare a 
balanced fermionic spin mixture with total atom numbers 
between 1.5 x 10 5 and 2.0 x 10 5 , with 10% systematic un¬ 
certainty. A mp = —9/2, —5/2 (mp = —9/2, —7/2) mix¬ 
ture with temperatures of 16(2)% (13(2)%) of the Fermi- 
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FIG. 1. Noise correlations (a) Schematic view of the ionic 
Hubbard model on a honeycomb lattice at half-filling. Circles 
denote lattice sites A and B, where larger circles indicate 
lower potential energy. The phase diagram exhibits two lim¬ 
iting cases: For A > [/,f a CDW ordered state is expected 
with two fermions of opposite spin (red, blue) on lattice sites 
B, and empty sites A. In the other limit (U A, t) a MI with 
one fermion on each lattice site should appear, (b) Measured 
noise correlation pictures obtained from absorption images of 
the atomic momentum distribution. Comparing panel 1 with 
panel 2, additional correlations appear due to broken inver¬ 
sion symmetry in the CDW ordered phase. When introducing 
strong interactions, these correlations are not observed any¬ 
more (panel 3), and the inversion symmetry of the density 
distribution does no longer reflect the broken inversion sym¬ 
metry of the lattice potential. Below each panel horizontal 
and diagonal cuts of the noise correlation image are shown. 
For the three different ratios of A and U , between 165 and 
201 measurements were taken each. We show the average of 
C(d x , dy) and C(d x , —d y ), which reflects the symmetry of the 
system. 


temperature, is then loaded into a three-dimensional 
optical lattice within 200 ms. Using interfering laser 
beams at a wavelength A = 1064 nm we create a hon¬ 
eycomb potential in the xy-plane, which is replicated 
along the z-axis [15, 31]. All tunnelling bonds are set to 
t/h = 174(12) Hz. The tunable lattice allows us to inde¬ 
pendently adjust the energy-offset A = [0.00(4),41 (1)]t 
between the A and B sub-lattice [32]. Depending on 
the desired interaction strength we either use the Fesh- 
bach resonance of the mp = —9/2,—7/2 mixture or the 
mp — —9/2, —5/2 mixture. 

We probe the spatial periodicity of the density distri¬ 
bution in the interacting many-body state by measuring 
correlations in the momentum distribution obtained after 


time-of-flight expansion and absorption imaging [33-38]. 
After preparing the system in a shallow honeycomb lat¬ 
tice with a given U and A, we rapidly convert the lattice 
geometry to a deep simple cubic lattice. This ensures 
that we probe correlations of the underlying density order 
rather than a specific lattice structure. The atoms are re¬ 
leased from the lattice and left to expand ballistically for 
10 ms. We then measure the density distribution, which 
is proportional to the momentum distribution of the ini¬ 
tial state n(q). From this, we compute the correlator of 
the fluctuations of the momentum distribution [33-39] , 

r( ,, = /(n(q 0 - d/2) ■ n(q 0 + d/2))dq 0 _ 

[> /(n(qo-d/2))(n(qo+d/2))dq 0 ’ U 

where the () brackets denote the statistical averaging over 
absorption images taken under the same experimental 
conditions. 

Owing to the fermionic nature of the particles, this 
quantity exhibits minima when d = m27r/A, with m a 
vector of integers [32]. This is illustrated by the anti¬ 
correlations of a repulsively interacting, metallic state 
with U = 4.85(9 )t and A = 0.00(4)£, shown in Fig. lb, 
left panel. There, the spatial periodicity of the atomic 
density follows the structure of the lattice potential, and 
minima in the correlator are observed for m = (0, ±2) 
and m = (±2,0). For A = 39.8(9)£, additional minima 
are observed at m = (±1, ±1), see Fig. lb, central panel. 
For a simple cubic lattice potential of periodicity A/2, the 
amplitude of these minima is given by [32] 


( 2p 27tA ( N a -N b ) 

V A’ \) a (7 V a ±A b ) 2 ' 


Thus, the observation of additional minima confirms the 
presence of CDW-ordering with Na ^ -/Vb- Finally, for 
A = 20.3(5)£ and U = 25.3(5)£, these additional minima 
are not observed any more (see Fig. lb, right panel), sig¬ 
nalling that with repulsive on-site interactions, the den¬ 
sity distribution does not reflect the externally broken 
inversion symmetry. In this case the interactions sup¬ 
press the CDW-order, despite the presence of a large A. 

Based on these measurements we expect the local dis¬ 
tribution of atoms on each lattice site to depend on 
the exact values of U and A. We measure the fraction 
of atoms on doubly occupied sites D using interaction- 
dependent rf-spectroscopy [9] . The number of doubly oc¬ 
cupied sites compared to the number of singly occupied 
sites is directly related to the nature of the insulating 
states [40, 41]: the MI state is signaled by a suppressed 
double occupancy while the CDW order is formed by 
atoms on alternating doubly occupied sites. 

In the experiment we set an energy-offset A and mea¬ 
sure D for different attractive and repulsive interactions 
U = [—24.6(13), ±29.1(7)]£. Fig. 2a shows D as a func¬ 
tion of U at constant A = 16.3(4)£. For strong attractive 
interactions we observe a large fraction of doubly occu¬ 
pied sites, which continuously decreases as U is increased. 
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FIG. 2. Double occupancy measurement (a) The measured 
double occupancy D as a function of the on-site interaction 
U for a fixed energy-offset A =s 16.3(4)t. (b) For different 
values of A (different colors) we obtain the double occupancy 
for a range of interactions U — [—24.6(13), 29.l(7)]t. Hol¬ 
low (full) circles represent attractive (repulsive) interactions. 
Vertical error bars show the standard deviation of 5 measure¬ 
ments and horizontal error bars the uncertainty on our lattice 
parameters. 

When tuning from attractive to weak repulsive interac¬ 
tions (A £/), we still observe a large D as expected for 
the CDW. For strong repulsive interactions (U A) the 
measured double occupancy vanishes, the density pat¬ 
tern no longer reflects the broken inversion symmetry of 
the lattice, confirming the suppression of the CDW or¬ 
dering. Fig. 2b shows D as a function of the energy scale 
V— A, which is the energy difference of a doubly occupied 
site neighbouring an empty site compared to two singly 
occupied sites in the atomic limit. For the largest nega¬ 
tive value of U — A we observe the highest D for all A. 
For positive values of V — A the double occupancy con¬ 
tinuously decreases and vanishes for the largest positive 
U — A, consistent with a MI state. In contrast, for the 
intermediate regime the measured D depends on the in¬ 
dividual values of U and A, as now the finite temperature 
and chemical potential itself play an important role and 
a detailed analysis would be required for a quantitative 
understanding, however we can qualitatively compare the 
dependence of D to an atomic limit calculation [32]. 

A characteristic feature of the MI and band insulating 
CDW state is a gapped excitation spectrum, which we 
probe using amplitude-modulation spectroscopy [9, 42]. 
We sinusoidally modulate the intensity of the lattice 
beam in y-direction by ±10% for 40 ms. Since the hon¬ 
eycomb lattice is created from several beams interfering 
in the xy-plane [15], this leads to a modulation in tunnel 
coupling t y of 20% and t x of 8%, as well as a modulation 
of U by 4% and A by up to 6%. The interlayer tun¬ 
nelling t z is not affected meaning that excitations only 
occur in the honeycomb plane. We set U = 24.4(5)£ and 
measure D after the modulation for frequencies up to 
v = 11.6 kHz (« 67 1). All measurements are performed 
in the quadratic-response regime [43] . 

Fig. 3a shows the measured spectra for different values 



Frequency hvlU A /U 

FIG. 3. Modulation spectroscopy measurement (a) Excita¬ 
tion spectra observed by measuring the double occupancy D 
from amplitude modulation spectroscopy of the lattice beam 
in y-direction for different energy-offsets A at repulsive on¬ 
site interaction U = 24.4(5)t. Solid lines are multiple Gaus¬ 
sian fits to the modulation spectra, (b) Schematics for the 
relevant energy scales \U — A| and U + A as a response to 
the lattice modulation, (c) Modulation spectroscopy of the 
lattice beam in z-direction. The measured excitation frequen¬ 
cies are shown as a function of A and compared to the value 
of U — 24.4(5)t (horizontal line). The inset shows the spa¬ 
tially dependent excitation spectrum, (d) Comparison of the 
measured excitation resonances (points) with the values of 
\U — A|, U ± A (lines). The area of the marker indicates 
the strength of the response (peak height) to the lattice mod¬ 
ulation. Full (empty) circles represent a positive (negative) 
response in double occupancy. Error bars as in Fig. 2, vertical 
error bars in (c),(d) show the fit error for the peak position. 


of A. The MI state exhibits a gapped excitation spec¬ 
trum, which is directly related to a particle-hole excita¬ 
tion with a gap of size U [9, 31, 43]. In the limit of A = 0 
we detect this gap as a peak in the excitation spectrum 
at v = U/h. With increasing A the single excitation 
peak splits into two peaks corresponding to different ex¬ 
citation energies [44]. The nature of the excitations can 
be understood as follows: The transfer of one particle 
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FIG. 4. Theoretical result for the kinetic energy response 
function xiy) °f the double occupancy on a modulated four 
site model as a function of A at constant U = 25 1 . Circu¬ 
lar (diamond) data points represent the response for the half 
filled (quarter filled) case. The area of the marker shows the 
relative size of the calculated response, whereas full (empty) 
data points have a positive (negative) response signal. 

costs approximately an energy of V — A if a double oc¬ 
cupancy is created on a B site and U ± A if it is created 
on an A site (see Fig. 3b). The excitation of additional 
double occupancies shows that atoms were initially pop¬ 
ulating both sub-lattices, as expected in the MI regime. 
For small A/U the system shows a clearly identifiable 
charge-gap, which vanishes if V ^ A. For large A the 
charge gap reappears, and a minimum in the spectra re¬ 
veals the breaking of double occupancies as a response to 
amplitude modulation. This is in agreement with the ex¬ 
pected band insulating CDW, where double occupancies 
are on the B sub-lattice and A sites are empty. 

The situation changes for amplitude modulation of the 
z lattice beam intensity by ±10%. In this case excitations 
are created along links perpendicular to the honeycomb 
plane. Since the honeycomb lattice is replicated along the 
z-axis, we observe a single peak at v = U/h, independent 
of the energy-offset A (see Fig. 3c). The inset of Fig. 
3c shows the direction dependent modulation spectrum 
for A = 8.5(2)t, which allows us to independently deter¬ 
mine the energy scales of the system in different spatial 
directions. 

We extract the excitation energies by fitting multiple 
Gaussian curves to our experimental data and compare 
our results with the values of \U — A|, U ± A and U in 
Fig. 3d. We observe a vanishing peak at U ± A for the 
largest A. This is expected as there are fewer and fewer 
atoms on A sub-lattice in the system for an increasing 
energy-offset. Our measurements are in good agreement 
with a picture based on nearest-neighbour dynamics. 

However, we observe additional peaks at v « U/h if 
U ~ A, which can not be understood in a two-site model. 
To rule out any higher-order contribution, we verified 
that the response signal has a quadratic dependence on 
the modulation parameters, as expected for quadratic re¬ 


sponse [43] . This additional peak was also observed in a 
purely 2D ionic Hubbard model (t z = h x 2 Hz), thus 
ruling out a contribution of excitations along the third 
direction [32]. 

To interpret the nature of the response at hv « U we 
calculate the kinetic energy response function 

X{v) = ^{m\5D\m)\(rn\K\0)\ 2 5(hv - e m0 ), (4) 

rn 

where the sum runs over all many-body states m, SD = 
D — (0|D|0) is the induced change in double occupancy, 
K = a an d e m o denotes the excitation energy 

measured above the ground state |0). We evaluate xiy) i n 
exact diagonalization of a cluster of four sites for varying 
filling fractions. 

The result shown in Fig. 4 for U/t = 25 clearly indi¬ 
cates that the peak at hv « U around U = A originates 
from regions of the lattice where the filling deviates from 
one particle per site [45]. In particular, for a configu¬ 
ration with two particles on four sites, the ground state 
at U = A is a configuration with negligible double oc¬ 
cupancy and only the lower sub-lattice sites are filled. 
The lattice modulation at hv ~ U then moves one par¬ 
ticle to an energetically costly site. For U = A, this 
configuration is resonantly coupled to a state where both 
particles are on the same, low-energy site. Hence, this 
process leads to an increase in the measured double oc¬ 
cupancy. The analysis of such a four-site cluster qualita¬ 
tively agrees with the observed signal at energy U in the 
intermediate (U ~ A) regime. 

In conclusion, we have realized and studied the ionic 
Hubbard model with ultracold fermions in an optical hon¬ 
eycomb lattice. Our observations show that increasing 
interactions suppress the CDW order and restore inver¬ 
sion symmetry of the density distribution. Additionally, 
we probed correlations beyond nearest-neighbour, which 
had not been accessible so far [46]. Future work can ad¬ 
dress open questions concerning the nature of the inter¬ 
mediate regime between the two insulating phases, which 
is theoretically debated and should depend on the dimen¬ 
sionality of the system [26, 47]. Furthermore, we can ex¬ 
tend our studies of the ionic Hubbard model to include 
topological phases by introducing complex next nearest 
neighbour tunnelling [30, 48, 49]. 
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SUPPLEMENTAL MATERIAL 
Preparation and optical honeycomb lattice 

To simulate the ionic Hubbard model we create a quan¬ 
tum degenerate cloud of 40 K by using a balanced spin 
mixture of the |F, mp) = |9/2, —9/2) and |9/2, —7/2) 
Zeeman states, which is evaporatively cooled in an op¬ 
tical dipole trap. Depending on the desired interac¬ 
tion strength we either use the Feshbach resonance of 
the mp = —9/2,—7/2 mixture to access an interaction 
range of U = [—24.6(13), 4.91 (9)] t or a mp = —9/2,—5/2 
mixture to reach strongly repulsive interaction strengths 
U = [11.7(2), 29.1(7)]£. The fermions, with temperatures 
between 16(2)% and 13(2)% of the Fermi-temperature, 
are then loaded into a three-dimensional optical lattice 
within 200 ms. The honeycomb lattice with staggered 
energy-offset A is created by interfering laser beams at 
a wavelength A = 1064 nm that give rise to the following 
potential [15]: 

V(x , y , z) = —cos 2 (kx + 6/2) — V x cos 2 (kx) 

—Vy cos 2 (ky) — V% cos 2 (kz) 

—2ayJVxVy cos (kx) cos (ky) cos </?, (5) 

where V ^ x yz are the s i n §l e beam lattice depths in 
each spatial direction, k = 27r/A, and the visibility 
a = 0.90(5). We interferometrically stabilize the phase 
cp = 0.00(3)7r and control the energy-offset A between 
the A and B sub-lattice by varying the value of 6 around 
7r. This is achieved by changing the frequency detuning 
between the X and the X (which has the same frequency 
as Y) beam. In the case of A = 0 the lattice depths are 
set to V x X Y z = [14.0(4), 0.79(2), 6.45(20), 7.0(2)]^ 
to prepare isotropic tunnelling bonds in the honeycomb 
lattice (t/h = 174(12) Hz), where Er is the recoil en¬ 
ergy. When breaking inversion symmetry (A % 0) we 
adjust the final lattice depths in order to keep t on all 
lattice bonds constant. Owing to the harmonic con¬ 
finement of the lattice beams and the remaining dipole 
trap our system has harmonic trapping frequencies of 
v x ,y,z = [55.6(7), 106(1), 57(1)] Hz. The lattice depths are 
independently calibrated using Raman-Nath diffraction 
on a 87 Rb Bose-Einstein condensate. For the measure¬ 
ment of the double occupancy D, both the independently 
determined offset in D of 2.2(3)% due to an imperfect 
initial spin mixture as well as the calibrated detection ef¬ 
ficiency of 89(2)% for double occupancies are taken into 
account [46]. 

Noise correlations 

The anti-correlations in the fluctuations of the momen¬ 
tum distribution stem from the fermionic nature of the 
particles. To illustrate this, let us consider the simple 


case of two identical fermions, occupying the lowest band 
of a one-dimensional optical lattice of periodicity A/2. 
Their wavefunctions can be expressed as Bloch waves, 
with a different quasi-momentum q for each particle, such 
that they obey the Pauli principle. If the momentum of 
the first particle is measured to be gi, then the second 
cannot be detected in any of the momenta q\ — nMn/X 
with n an integer. 

In the experiment, the momentum distribution 
n(g, r = 0) is accessed by suddenly releasing the atomic 
cloud from its confinement. After an expansion time suf¬ 
ficient to neglect the initial size of the cloud, the mo¬ 
mentum distribution has been converted to the spatial 
atomic density n(x,r) which is then directly measured 
by taking an absorption image. These two quantities are 
related by 

n(x , r) = n{q = t = 0) W 

In the following we use n to designate the atomic density 
in momentum space. In general, we are interested in the 
probability to detect two particles at momenta q\ and q^ 
simultaneously, which is given by 

P(<h,Q 2 ) = (' n(qi)n(q 2 )) - (n( qi )) (n(q 2 )) 

a,/3 

~ (^(«l)^a(</l)>(^(®)^(«2)), (7) 

where 4^ (g) creates and ^ a (q ) annihilates a particle 
with internal state a at momentum g. Representing the 
operators by their underlying Bloch-wave structure we 
have 

P(q 1 ,q 2 ) = \W(q 1 )\ 2 \W(q 2 )\ 2 Y J E e iA / 2 (fli(*-"*)+®(»-")) 

a,(3 k,l,m,n 

X bpj t>a,m b,3,n) — (^,fc b a ,m)(bp t i fys.n)) 

+ %i - 92) (8) 

a 

where k (b a ,k) creates (annihilates) a particle at the 
site k with an internal state a. Here, W(q) is the slowly 
varying envelope of the Bloch wave. The last term con¬ 
cerns the correlation of an atom with itself, which we are 
not considering here. To calculate the four-operator ex¬ 
pectation value, we assume that the atomic distribution 
is well described by Fock states ((b^ k bpj} = Sk,i$a,prik)i 
where rij . c is the number of particles on site k. Thus, we 
have 

P( qil q 2 ) = -\W( gi )\ 2 \W(q 2 )\ 2 

x EE n ^ n «P A/2(ft_i)(91_92) ( 9 ) 

a k,l 

P(qo,d ) = -\W( qo + d/2)\ 2 \W(q 0 -d/2)\ 2 

x E E A/2(fc_i)d , (10) 

a k,l 



7 



FIG. 5. Modulation spectroscopy measurement for the two- 
dimensional ionic honeycomb potential. Excitation spectra 
observed by measuring the double occupancy D after si¬ 
nusoidal modulation of the lattice depth V y for an energy- 
offset A = 24.3(6)£ at constant repulsive on-site interac¬ 
tion U = 24.1(4)t, thereby realizing the intermediate regime 
\U — A| ~ t. Although the tunnelling perpendicular to the 
honeycomb planes is tuned below 2 Hz we still observe the 
response at energy « U. Error bars show the standard devi¬ 
ation of at least 3 measurements. 


where we have introduced the center of mass go = (gi ± 
g2)/2 and the relative position d = q\ — The slowly 
varying dependence in go can be rewritten in terms of the 
momentum distribution 


(n(qi))(n(q 2 )) = \W(q 0 + d/2)\ 2 \W(q 0 - d/2)\ 2 N 2 (11) 


with N the total atom number. Thus, the correlations in 
momentum are fully characterized by 

r(H) = _ Jdq 0 P(qo,d) _ 

/ dq 0 (n(q 0 + d/2))(n(q 0 - d/2)) 

_ _ n a,k n a,l A/2 (k-l)d ^2) 

OL k,l 

We now show that this quantity is not only sensitive 
to the periodicity imposed by the lattice, but can also 
reveal underlying order in the density distribution. Con¬ 
sider the case where, for each internal state, the density 
takes the value nj^/M on even-numbered sites and n^/M 
on odd numbered sites, with M the number of internal 
states. The correlation signal then takes on the form 


C(d) 


n\ ± ± 2 n a n b cos(A/2 d ) 

M N 2 



A (k-l)d 


(13) 

The sum is equal to N 2 if d = m 2ir/X with m an integer 
and zero otherwise, and the correlation signal is then 
simply 


C(m 2 tt/X) 


(n a + (—l) m n B ) 2 

M 


(14) 


Thus, anti-correlations always appear at momenta 2m x 
27r/A, corresponding to the reciprocal lattice vector, 


while anti-correlations at momenta (2ra±l) x 2tt/X signal 
a staggering of the atomic density between the even- and 
odd-numbered sites. While this derivation was carried 
out for a one-dimensional lattice, it can readily be gener¬ 
alized to higher dimensions, provided the full information 
on the momentum density can be accessed. 


In the experiment, we measure the momentum distri¬ 
bution of the absorption images following a three step 
detection protocol. After preparing the system in a shal¬ 
low honeycomb lattice with a given U and A, the lattice 
depth is suddenly increased in 1 ms, which prevents any 
further evolution of the atomic density distribution. Sub¬ 
sequently, the lattice geometry is converted to a simple 
cubic lattice within 1 ms. Measuring the density distri¬ 
bution in the honeycomb lattice would lead to additional 
peaks at m = (±1,±1) due to the displacement of the 
lattice sites with respect to a square lattice. We can 
estimate the strength of these additional peaks with a 
simple model for a hexagonal lattice with A = 0 by plac¬ 
ing Gaussian wave packets at the position of each lat¬ 
tice site of the real potential. By calculating the Fourier 
transform of this system we find that the strength of the 
m = (±1,±1) peaks is a factor of 6 smaller than the 
minima of the correlator at position m = (0, ±2) and 
m = (±2,0). Therefore ramping to a simple cubic con¬ 
figuration ensures our observable probes correlations of 
the underlying density order rather than a specific lat¬ 
tice structure. Finally, the strength of the interactions is 
reduced within 50 ms using the Feshbach resonance, the 
atoms are released from the lattice and left to expand 
ballistically for 10 ms after which we measure the density 
distribution by absorption imaging. 


From this measurement, we compute the quantity C, 
which is shown in Fig. 1. As our imaging technique in¬ 
tegrates the density along the line of sight, we do not 
have access to the full information, but rather to the 
column density. Thus, the derivation presented above 
should be generalized to two dimensions, while the oc¬ 
cupancy along the third direction can be treated as an 
internal degree of freedom. Accordingly, the depth of the 
anti-correlation minima for a two-component Fermi gas 
will be divided by 2 N z , where N z is the typical num¬ 
ber of sites populated along the integrated direction. To 
achieve optimal signal-to-noise, we only consider atomic 
densities above 20 % of the maximum density. Further¬ 
more, we remove short-range correlations induced by the 
readout noise of the CCD chip of the camera by convo- 
luting the density distribution with a Gaussian of width 
A q = fc/25. Finally, we take advantage of the reflection 
symmetry of the momentum distribution, and average 
together C(d x ,d y ) with C(d x ,—d y ). Note that, by defi¬ 
nition, C(d x ,dy) = C(—d x1 —d y ). 
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FIG. 6. Theory comparison of the density per site in the 
trapped system and calculation of the double occupancy. For 
the calculation we use an entropy per particle of 1.5 ks (as 
measured before loading the atomic gas into the lattice) and a 
total atom number of 190000. (a) Density per site calculated 
for U = 25 1 and A = 0 using a local density approximation 
to include the harmonic trap. The high-temperature series 
(HTS) expansion up to second order of the grand canonical 
partition function is shown in grey. The atomic limit cal¬ 
culation (high-temperature series expansion in 0th order) is 
plotted in blue, (b) For different values of A (different colors) 
we calculate the double occupancy in the atomic limit using 
a local density approximation. We can compare the atomic 
limit calculation (blue) with a high-temperature series expan¬ 
sion up to second order (grey) for A = 0. Points show the 
measured double occupancy as plotted in Fig. 2b. 


Modulation spectroscopy 

The tunnelling t z between sites of adjacent layers can 
be controlled via the lattice depth Vjp To realize the 
ionic Hubbard model in two dimensions we suppress the 
tunnelling t z below 2 Hz by setting the lattice depth 
along the z-direction to = 30(1)£% thereby decou¬ 
pling the honeycomb planes. For the modulation spec¬ 
troscopy measurement we follow the same procedure as 
described in the main text and sinusoidally modulate the 
amplitude of the lattice depth in y-direction by ±10%. 
As a result we exclude possible contributions to the re¬ 
sponse signal, which may result from a residual coupling 
to the orthogonal direction. This ensures that the linear 
response measurement only probes energy scales that are 
realized within the xy-honey comb plane. Fig. 5 shows 
the excitation spectrum for the two-dimensional case. 
Even with suppressed tunnelling t z we observe a clear 
peak for modulation frequencies hv = U . As a result we 
can exclude that our response signal is resulting from an 
imperfect orthogonality of the lattice beams. 


Atomic limit calculation and local density 
approximation 

We can qualitatively estimate the density per site, tem¬ 
perature and double occupancy of our system from an 
atomic limit calculation (high-temperature series expan¬ 
sion of the partition function in 0th order) including the 
harmonic trap via a local density approximation (i.e. lo¬ 
cally varying the chemical potential). While this is not 
an exact theory capturing all details and not suitable for 
a precise experiment-theory comparison, this calculation 
gives a qualitative estimate for the typical system param¬ 
eters. In the case of A = 0 we can directly compare the 
obtained results for the density of atoms per site with 
a high-temperature series expansion up to second order 
[31] of the grand canonical partition function, see Fig. 6a. 
As expected, when comparing the results for the density 
of atoms per site we see a small deviation between the 
two theories. 

Using the atomic limit calculation with the local den¬ 
sity approximation we can estimate the fraction of half¬ 
filling and quarter-filling in our trapped system. Approx¬ 
imately a fraction of 46% of the atoms are within 10% 
of half-filling for values U = 25£ and A = Ot. The na¬ 
ture of the system totally changes for U = A = 25£, as 
a large fraction of the atoms (45%) consists of a quarter 
filled region. Deep in the two limiting cases (U A and 
A >> [/), we expect the approximation to provide better 
results. Given an entropy per particle of 1.5 ks (as mea¬ 
sured before loading the atomic gas into the lattice), and 
assuming adiabatic loading into the lattice, we find for 
the calculated temperature T values 3.3£ (5.5£), for pa¬ 
rameters of U = 30 1 and A = 0 (U = lOt and A = 40 1). 
As A and U are much larger than £, this confirms that 
the temperature is well below the charge gap for the MI 
and CDW states. 

For a qualitative comparison of our measurements of 
the double occupancy in Fig.2 we can calculate the frac¬ 
tion of double occupied sites D for various values of U 
and A using our simplified model at constant entropy per 
particle of 1.5As can be seen in Fig. 6b our results 
qualitatively confirm our observations that for increasing 
U the double occupancy is reduced well before reaching 
the point U = A. This calculation shows that the major 
contribution of the shift is resulting from a change in the 
chemical potential and temperature for different values 
of U and A. For A = 0 we can estimate the importance 
of the tunnelling t by directly comparing the double oc¬ 
cupancy obtained from the atomic limit calculation with 
a high-temperature series expansion up to second order. 




